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Abstract 

We present numerical results that allow a precise determination of the transition point and of 
the critical exponents of the 41) Edwards- Anderson Spin Glass with binary quenched random 
couplings. We show that the low T phase undergoes Replica Symmetry Breaking. We obtain 
results on large lattices, up to a volume V = 10^: we use finite size scaling to show the relevance 
of our results in the infinite volume limit. 

1 Introduction 

The question of whether short range Edwards Anderson (EA) spin glasses share the remarkable 
features of the infinite range Sherrington Kirkpatrick (SK) model is still an open one. In this work 
we present Monte Carlo simulations of the 4D EA Ising Spin Glass @, ||, H, ||, 0, H with a bimodal 
distribution of the quenched couplings, performed on large lattice volumes (thanks to the tempering 
and parallel tempering simulation technique |^, 0]), with a large number of samples, and down to 
low values of the temperature T. In this way we are able to obtain detailed information about the 
nature of the transition, to determine with good precision critical temperature and exponents, and to 
give strong evidence supporting the fact that the low-temperature phase is mean-field-like. A great 
deal of effort has gone in ensuring reliability of the data on delicate issues such as thermalization 
^ checks and consistency of data analysis. 

The paper is organized as follows: first of all we describe the model and the parameters of our 
MC simulation. We then present data related to the Binder cumulant and to the determination of 
Tc and i^. By analyzing the overlap susceptibility we determine the value of i]. Finally we discuss in 
detail about the probability distribution of the overlap P{q)- We present, among others, evidence 
for non-triviality of single sample Pj{q) and for a non-zero value of the position of the maximum of 
P{q), Qmax, in the thermodynamic limit. 

*e-niail address:Enzo. Marinari@ca.infn.it 
'''e-mail address: Francesco.Zuliani@ca.infn.it 
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2 The Numerical Simulation 



It is very difficult to run reliable numerical simulation of finite dimensional spin glasses. The main 
reason for such difficulties is the presence of many meta-stable states (responsible for aging effects 
as well as for many other peculiarities of spin glasses The Monte Carlo dynamics gets easily 
trapped, and the system only probes a restricted part of phase space. 

Many algorithmic solutions have been proposed to improve the speed of thermalization of these 
systems. All these techniques are related to density scaling methods (see [|10| for a review and 



references): we use here the maybe simplest implementation of these ideas, the parallel tempering p], 
where a number of configurations of the system are allowed to exchange their temperature (for 
multi-canonical methods, that are strongly related and have in principle an even wider range of 
applicability, see for example []TT|). Thanks to parallel tempering we have been able to thermalize 
systems of volume V = 10'* down to T ^ 1.2(0.6Tc). 

We study the AD Edwards-Anderson Ising spin glass with binary couplings, with Hamiltonian 

H =-J2 Jij(^i(^j , (1) 

where the sum runs over nearest neighboring sites, the cTj are ±1 Ising spins, and the couplings are 
quenched variables drawn with probabihty | among the two values { — 1, +1}. The overlap among 
two different systems is defined as 

?"'^ = ^E<^f> (2) 

i 

where a and (3 denote two configurations of the system in the same realization of the quenched 
disorder. The overlap probability distribution for a given sample is 

P,(9) = (%-g"'^)) , (3) 
where (...) denotes the usual Gibbs average. Its average over samples is 



P{q) = Pj{q) , (4) 

and its moments are defined as 



gin) 



{qn) = Jdqq- P{q) . (5) 



We always denote by (■ ■ ■) the thermal averages and by the disorder averages. 

Our simulation have been performed on a set of workstations, using a multi-spin-coding program 
that was inspired by the work of ||12|. We have selected the parameters of our Monte Carlo and 



parallel tempering runs such to guarantee a complete thermalization of the measured observables. 
We will discuss this issue in some detail. 

Table (0) summarizes the relevant parameters used in the simulation: we give among others the 
number of thermalization steps, of measurements steps, the number of different disorder realizations 
and the temperature ranges investigated by tempering. Temperature values have been chosen 
uniformly spaced in the interval between Tmin and Tmax- 

We have used different methods to verify that we have correctly thermalized the systems. Using 
"parallel tempering", one is actually performing a generalized Markov chain where systems at 
different temperatures are allowed to "move" in temperature-space too. A necessary condition for 
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Table 1: Parameters of the Tempered Monte Carlo runs. 

the Markov chain to be effective in de-correlating different measurements is the fact that each system 
spans at least a few times all the allowed temperature range during the simulation. In this respect 
we check a posteriori that the probability of swapping temperature has been of order 0.5 (ensuring 
in this way that a single system did not get stuck at a specific value of T) and that the histogram 
counting the time that each system has spent at each temperature is fairly flat. This requirement 
is fulfilled in all our simulations, for all T and L values. 

Another very strong check of thermalization is the fact that the single sample Pj{q) are sym- 
metric in the limits of the statistical significance of the histogram. This is very well verified as can 
be seen for example in figure O where we plot Pj{q) for selected samples. 



3 The Binder Parameter, Tc and u 

We start by discussing the overlap Binder parameter. We will use it to qualify the phase transition, 
and to determine the critical temperature and the first of the critical exponents, v. We will use and 
describe different methods to compute the quantities we are interested in. Our statistical sample of 
configurations is a large sample, and our set of data precise (even as far as the dependence over the 
lattice volume V is concerned): we will show that different analysis styles give compatible (precise) 
results. 

We define the usual overlap Binder parameter as 
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(6) 



The Binder parameter is an adimensional quantity, and its value at the critical point is universal. 
Close to Tc its leading behavior is 

g{L,T)c^g{L^iT-T,)) . (7) 

In usual ferromagnetic systems the infinite volume limit of the magnetization Binder cumulant is 
in the warm phase (where the distribution of the order parameter is Gaussian) and 1 in the broken 
phase: for a spin glass with replica symmetry breaking and hence a non-trivial distribution of the 
overlap order parameter, the transition is signaled by a non-trivial value of g in the broken phase 
(in the warm phase one expects an infinite volume limit of zero). In both cases the location of 
Tc is signaled by the crossing of the curves of g versus T for different values of the lattice size L 
(asymptotically for large L): large L curves are lower for T > and higher for T < T^. We show in 
figure |I| g versus T for different L values. The crossing point is close to T ^ 2 for all lattice values. 
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Figure 1: Binder parameter g versus T, for different values of the linear lattice size L (see caption 
in the plot). 



and the value of the Binder cumulant at criticality, gc is close to 0.45. Also error analysis has been 
a sensitive issues. We have always used a jackknife or a bootstrap error analysis [|1^] directly on 
the fitted parameters to determine errors. Still one has to keep in mind that statistical errors come 
together with systematic errors, due to the functional form one decides to try to fit (typically the 
asymptotic scaling form, that on finite size lattices is affected by power corrections). The two types 
of errors have to be kept under control separately. 

In figure |l| the crossing of the different gL curves is very clear. It is interesting to stress the 



difference with the three dimensional case |T^ , where the crossing at looks more like a merging 



of the different curves. 3d is (very) close to the lower critical dimension, while in Ad we are in a safe 
region: potentially this is important to make the physical picture easier to understand. 

Let us discuss a first naive approach to the data. By looking at the crossing of the curves giiT) 
versus T for different (L, L + 1) values one sees that one cannot extract a systematic dependence of 
the crossing point (and hence of the estimate of the effective critical temperature Tc(L, i^+1)) over L. 
Any systematic trend is smaller than the statistical error (maybe just showing a systematic average 
decrease of the estimate Tc{L, L + 1) when going from smaller to larger L values). The preferred 
value of Tc is slightly larger than 2.00. A first naive estimate of u can be done by linearizing 
gL{T) around the estimate we have given for and by evaluating the logarithm of the slope ratio 
(divided by the logarithm of the two lattice sizes ratio, \og(^j^^). With this method, one gets 
a first estimate for a set of effective exponents z^(L,L + 1). Here too, one cannot distinguish any 
clear strong dependence over the lattice size: the error one gets on u is completely correlated to 
the variation of the estimate of Tc. For larger Tc one estimates a lower value of u, while for lower 
estimates of Tc one gets larger estimates for ly. The error is dominated by this effect. The estimate 
for u is close to 1. 
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To get a reliable estimate of and of v we have used two methods of analysis of g (see for 
example the discussion of the analysis of reference [0, In the first approach we linearize the 
data close to (for all L values) and we run a global fit to all data: we fit and v for the 
two variable function giiT) (as we said, linearized close to Tc). We use data in a T range around 
the interval 1.9 — 2.1. We estimate the errors over the fit parameters (Tc and u) by a jackknife 



approach [|T3[: we repeat the fit approximately K times over a subsample of the data containing all 
of our statistical sample but a fraction The error is estimated by looking at fluctuations of the 
results of the K fits, and by accounting for the fact they are correlated [|l^. We also repeat the 



fits by discarding the smaller L values, to check if we can observe any systematic drift (again with 
good accuracy the average value of the result does not seem to depend systematically over the L 
range selected). Results are very stable, and the value we estimate for u systematically comes out 
to be close to 1.10. 

In the second approach, that comes in different flavors, one only uses data in the warm phase. 
This method leads to a smaller statistical error, that is balanced from a larger systematic incertitude 
(since we only select data at a given distance from Tc, and approaching Tc leads to a systematic 
drift of the estimate). In this case we start by selecting a threshold value for g, g* < gc- We start 
with low values of g*, and we approach gc from below: we cannot get too close to gc or the merging 
of the curves for different sizes makes the error over the measurement too large (we use values of 
g* going from = 0.2 to 0.4. We use a polynomial fit to interpolate the data for g{T), at different L 
values. We have decided to use a polynomial of degrees four (we have checked it guarantees stable 
fits and consistent results), and we fit a T range in the critical region (for L = 3 we use the data 
in the T range 1.5 — 2.8, for L = 10 we use the range 1.88 — 2.16). We define now Tc{L, g*) as the 
crossing point of the fitted polynomial with the horizontal line at g*, and v* as 

limrc(L,/)=Tc(/) + -4. (8) 

When g* ^ gc J^* ^ ^- If g* is too small violations of scaling are dominant, while if one approaches 
too much gc the merging of the g curves makes the error over the determination of v* overwhelming. 
The errors have been estimated by using a bootstrap approach (very similar in spirit to the jackknife 
technique, see |1^): one emulates fake sets of data with a Gaussian distribution around the real 
measurements, fits these multiple sets of fake data and compute the errors over the fit parameter. 



We note at last that we have also used a variation of this second method, described in |n 
based on the direct analysis of the derivative of g with respect to T. Also this method gives results 
that are compatible with the other ones. 

Our final estimates, averaged over the results obtained using these different approaches, are 

Tc = 2.03 ± 0.03 , (9) 

and 

= 1.00 ±0.10 . (10) 

In the rest of this paper we will use these two values as our best estimates of Tc and ly. We show 
in figure || the data for giiT) rescaled by using these two values: the scaling turns out to be very 
satisfactory. 
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Figure 2: versus Li(T — Tc), with z/ = 1.0, Tc = 2.03. 

4 The Overlap Susceptibility and 77 

The determination of the overlap susceptibihty, Xg, provides various possible ways to determine the 
exponent rj (and hence of the exponent 7) . In a spin glass in the RSB phase the overlap susceptibility 



is expected to diverge for all values of T <Tc. We show Xq versus T in figure ^. 
The first method is based on the fact that we expect that at T = Tc 



(12) 



We use a linear interpolation of the data in the region close to T^. As in the case of g the error 
in the estimate is mainly related to the choice of T^. The fit at T = 2.03 by using L > 3 gives an 
estimate of 0.28. 

In the second method we use data where L ^ ^. We go as close to Tc as possible, under the 
condition that data on our larger lattice (L = 10) coincide, in our statistical accuracy, with the ones 
at L = 8. Here we expect that 



X,(T)^(T-T, 



(13) 



We can use data down to T = 2.5 (i.e. at a AT ^ 0.5 from T^), where finite size effect start to be 
sizable even at L = 10. We show our best fit (in a T interval of = .2) in figure In this region 
we have a stable fit, with rj close to —0.4. Even if this second measurement is not very precise 
(we have to stay quite far from the critical region) it is interesting the fact that we get a coherent 
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Figure 3: Overlap susceptibility Xq versus T, for different L values. 



determination of r/, by using a completely different scaling region than in the former analysis (the 
new analysis also depends on the value of v we have determined by using g). 

The last approach we use for determining rj is based on the analysis of the scaling properties of 
the distribution probability P{q) of the overlap order parameter q in the region g ~ at T = Tc. 
We analyze the behavior of P{q) in the next section, but we discuss now the scaling of -P(O) of Tc 
in order to define our determination of r]. At T = we expect 

P{q ~ 0) ~ , (14) 

2+77 

i.e. in c? = 4 a scaling with L~ . We find a very good best fit (we do not include the L = 3 data), 
with an r] value close to —0.3. 

By considering all the methods we have discussed in this section we give our final estimate 

= -0.30 ± 0.05 , (15) 

that we will use in the rest of our analysis. In figure |^ we plot Xq rescaled by using our best fits. 
The rescaling works fine. 

Let us also notice that we have a good agreement with the results reported in 0] for the Ad 
EA model with Gaussian couplings. There the authors find p ~ 1.06, and t] ~ —0.35. Universality 
seems to work. 



5 P{q) 

In the two former sections we have shown that the AD EA model undergoes a phase transition, and 
we have determined its location and the critical exponents. Now we will try to qualify it in better 
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Figure 4: Best fit to the overlap susceptibility Xq versus T — at T > Tp. 
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Figure 5: Rescaled overlap susceptibility -^S^ versus L^{T — Tc), with Tc = 2.03, v 
T] = -0.30. 



1.0 and 



8 





-1 -0.5 0.5 1 

Figure 6: P{q) at T = 1.2 (broken phase), for different lattice sizes. 

detail, by determining and analyzing the probability distribution of the order parameter, P{q). 

In figure ^ we show our average P{q) (averaged over the different disorder realizations) at 
T = 1.2 < Tc- When increasing the lattice size the peak where P{q) is maximum shifts to lower q 
values: for showing that there is a phase transition to a phase with a non zero expectation value 
of q we have to show that the peak does not go to g = when L oo. The plateau of P{q) for 
g ~ does not lower when increasing L, as we will discuss better in the following. We remind the 
reader that in the RSB Parisi Mean Field scenario the P{q) is (in zero magnetic field) a non trivial 
function, that in the infinite volume limit is formed by a 5 function at g = qEA and by a regular 
part that extends down to g = 0. On the contrary if the broken phase has the same structure of 
the one of an ordered ferromagnet P{q) has to become asymptotically the sum of two 6 functions 
at ±qEA- 

For sake of comparison we show in figure |^ what happens in the warm phase, where the average 
P{q) shrinks to a Gaussian distribution around q = when L — > oo. 

In figure |^ we compare P{q) at different values of T on the same lattice size. From the single 
peaked shape at high T one gets a clear double peaked structure, with a clear plateau at low q, in 
the low T region. 

As we have said, in order to establish that we are having a real phase transition in the infinite 
volume limit we have to show that the value of g = qmax where P{q) is maximum does not go to 
zero. We start by plotting in figure ^ qmax versus L~^-^ (the exponent 1.3 comes from our best fit, 
see later). It is easy to see that an asymptotic value qmax = seems unplausible. 

Making this last statement more quantitative needs a more careful analysis. In order to do that 
we fit 

qmax{L) = qmax{00) + — , (16) 




Figure 7: P{q) at T = 2.2 (warm phase), for different lattice sizes. 




Figure 8: P{q) at L = 8, for all values of T. From single peak in g = to continuous part and 
double peak at large q for increasing T. 
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Figure 9: q^ax versus L 
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Figure 10: q^ax versus L and our two best fits. 
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Figure 11: P{q ^ 0) versus L. T = 1.2. 



both with Qmaxioo) = and by allowing for it a non zero value. In figure we plot the values of 
Qmax versus L and the results of the two best fits, one with a fitted value of Qmaxioo) and the second 
with fixed qmaxi^o) = 0. This second fit is clearly unsuitable, and it has a very high value of x^- In 
the first fit we get 



(00) = 0.548 ±0.006 



;i7) 



that is our best estimate for the position of the 6 function at qea in the infinite volume limit. 
We estimate a = 1.3 ± 0.1 (in the fit with a zero asymptotic value one finds the very small value 



a 



0.2). 



In figure |TT| we show the value of P{q) close to g = (averaged over a small q range, where P{q) 
is remarkably constant, in order to diminish statistical fiuctuations) as a function of L. One cannot 
observe any statistically significant decrease of this value for increasing large lattice volume (there 
is a small decrease only for small volumes). The most plausible implication of this evidence is that 
the system has many stable states, and that the cold T phase is characterized by Replica Symmetry 
Breaking (even if it has to be stressed that this evidence is not as strong as the one implied by the 
figure |10D. 

In figure |l^ we plot Pj{q) for selected samples, at T = 1.2, L = 10. One can see here that they 
are very complex distributions: such a pattern is typically related to the presence of many states 
(it has to be notice however that the small side peaks are not always there because of the presence 
of a real state). 

To be more quantitative we have measured the percentage of disorder configurations such that 
Pj{q) has 1, 2, 4 and 6 peaks versus L, and we plot it in figure |T3|. The number of configurations with 
a complex phase space {Pj{q) with many peaks) increases strongly with L. We use this evidence to 
rule out the picture of a modified droplet model, that has been discussed, among others, in [ITB] and 
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Figure 12: Pj{q) for selected samples. T = 1.2, L = 10. 




Figure 13: Percentage of disorder configurations such that Pj{q) has 1, 2, 4 and 6 peaks versus 
L. The number of configurations with a complex phase space {Pj{q) with many peaks) increases 
strongly with L. 
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Figure 14: (g^)^, {q^} and (g^) vs. T. 



in references therein. The picture of the modified droplet models imphes that for each realization 
of the quenched disorder there are (in the cold phase) only two ground states, but that the value 
of Qea (i-e. the support of the 6 function that constitutes the Pj{q)) depends on the sample. Here, 
on the contrary, the number of states for a given sample is strongly increasing with L (and with 
decreasing T). 



6 Sum Rules 



In this section we discuss another important feature of the broken phase of the 4D EA model. The 
starting point for this analysis can be for example the relation: 

W = IW)' + IW) . (18) 

This is one of a set of relations that are valid in the Mean Field Theory of Spin Glasses 
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The work contained in |T8| has established numerically that these relations are satisfied with good 
accuracy also in finite dimensional spin glasses. Following these findings a rigorous and theoretical 
analysis has improved our understanding of such set of sum rules [|19|, ^ ^ : they are strongly 
related to the ultrametric properties of the phase space. 

First of all we show evidence that the relation (|18D has a non-trivial content in the low- 
temperature phase (in the high T phase it is satisfied in the form = 0). Figure |14| shows that the 
values of the three quantities involved in (0) are significantly different from zero below Tc (we have 
already shown in better detail that the infinite volume of such quantities is non-zero). 



In figure |I5| we show the difference between the left hand side and the right hand side of (p!q). 
The two contributions cancel out with good accuracy (to 3 significant figures), and asymptotically 
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Figure 15: Left hand side minus right hand side of equation (p!8[) vs. T. 



for large lattice size the difference extrapolates to zero. 

Another possible way to visualize the result is plotting the ratio of the left hand side and the 
right hand side. As figure |16| shows, for T below T^. we get identically one, while for T ^ oo we get 
the value |, expected for a Gaussian P{q). 



We have also fitted the difference plotted in figure |T5], for various values of temperatures T < T^: 
in all cases a fit to an asymptotic zero value with power corrections works very well, and the exponent 
of the corrections is close to 3 for all T values. As an example we plot the data together with the 
best fit for T = 1.4 in figure |17. 



7 Conclusions 

In this note we have been able to give strong evidence for mean field behavior of the Ad Ising spin 
glass with binary couplings. Life in the Ad model is easier than in Sd, where even after a large 
number of intense numerical simulations the evidence for a phase transition is still slightly marginal 
(even if, at this point, convincing enough). In our case already the crossing of the Binder cumulants 
is the very clear signature of a typical phase transition (as opposed to the quasi-merging, quasi- 
Kosterlitz-Thouless behavior of the 3(i case). It is clear that 3d is very close to the lower critical 
dimension, and that there observing the effects of the physical critical point is dramatically difficult. 
Ad is on the safer side, and numerical simulations show that very clearly. 

We have been able to determine critical exponents precisely, and to enter in the large volume 
region with good accuracy. For example we have been able to show that the peak of P{q) is not 
going to g = for increasing lattice size, and (with a slightly worst accuracy and level of reliability) 
that the plateau at g ~ does not decrease with increasing lattice size. Also we remind the reader 
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Figure 16: Ratio of left hand side and right hand side of equation (|18|) versus T. 
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Figure 17: Left hand side minus right hand side of equation (ITB]) versus L, at T = 1.4, and and 
best fit to a zero constant value with a simple power correction. 
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that non-trivial sum rules are satisfied with very good accuracy. So, thinks look quite clear in the 
id case. 
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